Optimizing Urban Logistics with Network Analysis in RotterdamΒΆ

This project explores how network analysis and geospatial data can be used to optimize logistics operations in an urban environment. Using real geospatial datasets of distribution centers in the Netherlands and the road network of Rotterdam, we apply Python libraries such as OSMnx, NetworkX, GeoPandas, and Folium to solve two common logistics problems:

  1. Shortest Route Optimization between two warehouses
  2. Closest Facility Assignment for delivery points to their nearest logistics hub
  3. Service Area Analysis
  4. Facility Location Optimization

The goal is to simulate real-world supply chain scenarios using data-driven routing strategies and provide interactive visualizations to better understand the results.

The dataset used in this analysis can be found in the following link: https://data.4tu.nl/articles/dataset/Dutch_Distribution_Centres_2021_Geodata/19361018

InΒ [1]:
# Data manipulation
# ==============================================================================
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as px

pd.set_option("display.max_columns", None)

# Data visualisation
# ==============================================================================
import geopandas as gpd
import osmnx as ox
import networkx as nx
import folium


# Load the GeoPackage file and check its structure
# ==============================================================================
gdf = gpd.read_file('distributioncenters2021NL_repo.gpkg')
display(gdf.shape)
display(gdf.columns)
display(gdf.head())
/opt/anaconda3/lib/python3.12/site-packages/pyogrio/geopandas.py:265: UserWarning: More than one layer found in 'distributioncenters2021NL_repo.gpkg': 'buildings2021NL' (default), 'estates2021NL', 'corridor2021NL'. Specify layer parameter to avoid this warning.
  result = read_func(
(26951, 28)
Index(['osm_id', 'BAG_id', 'Constr_yea', 'Year_cat', 'Footprint', 'Size_cat',
       'Function', 'Function_n', 'Job_dens', 'Job_dens_n', 'RIN', 'Year',
       'Name', 'Locality', 'Gross_ha', 'Net_ha', 'Max_cat', 'Rail_cat',
       'Water_cat', 'Muni_nr', 'Municip', 'Prov_nr', 'Province', 'Corop_nr',
       'Corop', 'Corridor', 'Current', 'geometry'],
      dtype='object')
osm_id BAG_id Constr_yea Year_cat Footprint Size_cat Function Function_n Job_dens Job_dens_n RIN Year Name Locality Gross_ha Net_ha Max_cat Rail_cat Water_cat Muni_nr Municip Prov_nr Province Corop_nr Corop Corridor Current geometry
0 273364927 0441100000002816 1972 1970s 778.057432 S Trade, Import, Export 1 medium 3 1741 1971.0 Witte Paal Schagen 76.0000 43.0000 4 A E 441 Schagen 27 Noord-Holland 18 Kop Van Noord-holland 0 0 MULTIPOLYGON Z (((115780.809 532845.953 0, 115...
1 279126121 0971100000022139 1978 1970s 734.322479 S Trade, Import, Export 1 medium 3 203678 1960.0 Kerensheide Stein 19.8223 16.9625 3 D E 971 Stein 31 Limburg 39 Zuid-limburg 0 1 MULTIPOLYGON Z (((182791.839 331932.418 0, 182...
2 615938117 1931100000000715 2018 2010s and later 836.482827 S Trade, Import, Export 1 low 2 171034 2012.0 Zevender Schoonhoven 15.0000 10.7000 4 A E 1931 Krimpenerwaard 28 Zuid-Holland 28 Oost-zuid-holland 0 1 MULTIPOLYGON Z (((118809.492 440320.413 0, 118...
3 284074624 0302100000119189 2013 2010s and later 2048.445693 S XXL retail center 3 low 2 1203 1943.0 Kom Nunspeet Nunspeet 8.7132 8.3971 5 A E 302 Nunspeet 25 Gelderland 13 Veluwe 0 1 MULTIPOLYGON Z (((181723.773 487021.56 0, 1817...
4 280745712 0050100000352829 1999 1990s 1610.136749 S Trade, Import, Export 1 very high 5 213 1975.0 Trekkersveld I en II Zeewolde 160.1670 93.9994 5 A A 50 Zeewolde 24 Flevoland 40 Flevoland 0 1 MULTIPOLYGON Z (((163088.239 484983.775 0, 163...
InΒ [2]:
# Filter the 'gdf' dataset to include 'Transport and Logistic services' and the city of Rotterdam, since it will be the focus of this analysis
# ==============================================================================
rotterdam = gdf[
(gdf['Function'] == 'Transport and logistic services') &
(gdf['Municip'] == 'Rotterdam')
]

# Check new dataframe structure
# ==============================================================================
display(rotterdam.shape)
display(rotterdam.head())
display(rotterdam.tail())

# Define the geographical location to download road network
# ==============================================================================
place_name = 'Rotterdam, Netherlands'

# Download the road network
# ==============================================================================
G = ox.graph_from_place(place_name, network_type= 'drive')

# Plot the street network
# ==============================================================================
ox.plot_graph(G)
(323, 28)
osm_id BAG_id Constr_yea Year_cat Footprint Size_cat Function Function_n Job_dens Job_dens_n RIN Year Name Locality Gross_ha Net_ha Max_cat Rail_cat Water_cat Muni_nr Municip Prov_nr Province Corop_nr Corop Corridor Current geometry
33 271881098 0599100000653535 1999 1990s 784.190901 S Transport and logistic services 2 very high 5 170371 2010.0 Europoort Rotterdam 1700.0 1369.0 5 B A 599 Rotterdam 28 Zuid-Holland 29 Groot-rijnmond 1 1 MULTIPOLYGON Z (((72649.646 436919.835 0, 7266...
75 273554871 0599100000643345 1990 1990s 2362.159063 S Transport and logistic services 2 very high 5 170143 1964.0 Prins Alexander Rotterdam-noord 35.0 29.0 3 A E 599 Rotterdam 28 Zuid-Holland 29 Groot-rijnmond 1 1 MULTIPOLYGON Z (((97149.757 440577.129 0, 9718...
181 554042637 None 1994 1990s 4874.941239 M Transport and logistic services 2 no data 99 170380 1973.0 Eemhaven Rotterdam 358.0 292.0 4 B A 599 Rotterdam 28 Zuid-Holland 29 Groot-rijnmond 1 1 MULTIPOLYGON Z (((88339.665 433065.91 0, 88329...
228 273326689 0599100010036685 1988 1980s 3352.465527 M Transport and logistic services 2 no data 99 170373 1970.0 Botlek Europoort-oost Rotterdam 1200.0 1069.0 5 B A 599 Rotterdam 28 Zuid-Holland 29 Groot-rijnmond 1 1 MULTIPOLYGON Z (((77162.415 434307.927 0, 7718...
229 273326686 0599100000754059 2001 2000s 10457.050191 L Transport and logistic services 2 low 2 170376 1990.0 Distripark Botlek Rotterdam 120.0 108.0 4 B A 599 Rotterdam 28 Zuid-Holland 29 Groot-rijnmond 1 1 MULTIPOLYGON Z (((76763.925 432037.599 0, 7677...
osm_id BAG_id Constr_yea Year_cat Footprint Size_cat Function Function_n Job_dens Job_dens_n RIN Year Name Locality Gross_ha Net_ha Max_cat Rail_cat Water_cat Muni_nr Municip Prov_nr Province Corop_nr Corop Corridor Current geometry
26041 273554710 0599100000617012 1976 1970s 1782.039078 S Transport and logistic services 2 very low 1 170143 1964.0 Prins Alexander Rotterdam-noord 35.0 29.0 3 A E 599 Rotterdam 28 Zuid-Holland 29 Groot-rijnmond 1 1 MULTIPOLYGON Z (((96511.132 440308.944 0, 9654...
26279 863262550 0599100110003461 2020 2010s and later 33920.430577 XL Transport and logistic services 2 no data 99 170341 2008.0 Maasvlakte Rotterdam 2728.0 2297.0 5 B A 599 Rotterdam 28 Zuid-Holland 29 Groot-rijnmond 1 1 MULTIPOLYGON Z (((58428.859 439331.876 0, 5851...
26280 863262551 0599100110004080 2020 2010s and later 20116.635052 XL Transport and logistic services 2 no data 99 170341 2008.0 Maasvlakte Rotterdam 2728.0 2297.0 5 B A 599 Rotterdam 28 Zuid-Holland 29 Groot-rijnmond 1 1 MULTIPOLYGON Z (((58522.069 439072.086 0, 5864...
26445 271523822 0599100000659860 2000 2000s 5113.122114 M Transport and logistic services 2 very high 5 170164 1971.0 Hoogvliet Hoogvliet 28.0 19.8 3 A E 599 Rotterdam 28 Zuid-Holland 29 Groot-rijnmond 1 1 MULTIPOLYGON Z (((85042.696 432049.426 0, 8506...
26896 271871679 0599100000047089 1987 1980s 803.993757 S Transport and logistic services 2 very high 5 170371 2010.0 Europoort Rotterdam 1700.0 1369.0 5 B A 599 Rotterdam 28 Zuid-Holland 29 Groot-rijnmond 1 1 MULTIPOLYGON Z (((71684.971 437914.726 0, 7170...
No description has been provided for this image
Out[2]:
(<Figure size 800x800 with 1 Axes>, <Axes: >)
InΒ [3]:
from scipy.spatial import cKDTree

# Convert the graph to GeoDataFrames
# ==============================================================================
nodes, edges = ox.graph_to_gdfs(G)

# Reproject Transport and logistics buildings to match network CRS
# ==============================================================================
rotterdam = rotterdam.to_crs(nodes.crs)

# Build KDTree from network node coordinates
# ==============================================================================
node_coords = np.array([(point.y, point.x) for point in nodes.geometry])
tree = cKDTree(node_coords)

# Get coords of logistics points (using centroids of polygons)
# ==============================================================================
building_coords = np.array([(geom.centroid.y, geom.centroid.x) for geom in rotterdam.geometry])

# Find nearest node index for each index
# ==============================================================================
distances, indices = tree.query(building_coords, k=1)

# Add nearest node ID to the 'rotterdam' dataframe
# ==============================================================================
rotterdam['nearest_node'] = nodes.iloc[indices].index

display(rotterdam.head(3))
display(rotterdam.columns)
osm_id BAG_id Constr_yea Year_cat Footprint Size_cat Function Function_n Job_dens Job_dens_n RIN Year Name Locality Gross_ha Net_ha Max_cat Rail_cat Water_cat Muni_nr Municip Prov_nr Province Corop_nr Corop Corridor Current geometry nearest_node
33 271881098 0599100000653535 1999 1990s 784.190901 S Transport and logistic services 2 very high 5 170371 2010.0 Europoort Rotterdam 1700.0 1369.0 5 B A 599 Rotterdam 28 Zuid-Holland 29 Groot-rijnmond 1 1 MULTIPOLYGON Z (((4.19026 51.91467 0, 4.19049 ... 44519439
75 273554871 0599100000643345 1990 1990s 2362.159063 S Transport and logistic services 2 very high 5 170143 1964.0 Prins Alexander Rotterdam-noord 35.0 29.0 3 A E 599 Rotterdam 28 Zuid-Holland 29 Groot-rijnmond 1 1 MULTIPOLYGON Z (((4.5457 51.95062 0, 4.54624 5... 44550916
181 554042637 None 1994 1990s 4874.941239 M Transport and logistic services 2 no data 99 170380 1973.0 Eemhaven Rotterdam 358.0 292.0 4 B A 599 Rotterdam 28 Zuid-Holland 29 Groot-rijnmond 1 1 MULTIPOLYGON Z (((4.41902 51.88213 0, 4.41887 ... 1582380654
Index(['osm_id', 'BAG_id', 'Constr_yea', 'Year_cat', 'Footprint', 'Size_cat',
       'Function', 'Function_n', 'Job_dens', 'Job_dens_n', 'RIN', 'Year',
       'Name', 'Locality', 'Gross_ha', 'Net_ha', 'Max_cat', 'Rail_cat',
       'Water_cat', 'Muni_nr', 'Municip', 'Prov_nr', 'Province', 'Corop_nr',
       'Corop', 'Corridor', 'Current', 'geometry', 'nearest_node'],
      dtype='object')

Shortest Possible RouteΒΆ

In this section, we calculate the shortest path between two randomly selected logistics centers using the Dijkstra algorithm. This is a classic logistics optimization problem focused on minimizing travel distance across a real-world road network.

Key Objectives:

  • Learn how to compute optimal routes between fixed locations.
  • Understand how to use graph theory (via NetworkX) to traverse road networks.
  • Visualize the optimal route on an interactive map using folium.

This analysis is foundational for route planning, helping logistics companies reduce fuel costs and delivery times.

InΒ [4]:
# Define warehouse origin and end point by index (select two random warehouses)
# ==============================================================================
origin_node = rotterdam.iloc[0]['nearest_node']
destination_node = rotterdam.iloc[20]['nearest_node']

# Convert node IDs to lat/lon coordinates
# ==============================================================================
origin_coords = (G.nodes[origin_node]['y'], G.nodes[origin_node]['x'])
destination_coords = (G.nodes[destination_node]['y'], G.nodes[destination_node]['x'])

# Find the shortest path and the shortest path length using NetworkX (apply Dijkstra's algorithm)
# ==============================================================================
shortest_path = nx.shortest_path(G, origin_node, destination_node, weight= 'length')
path_length = nx.shortest_path_length(G, origin_node, destination_node, weight='length')

print(f" Shortest Path: {path_length:.2f} methers")

# Visualise using Folium 
# ==============================================================================
m = folium.Map(location=origin_coords, zoom_start=14)

folium.Marker(
    location=origin_coords, 
    icon=folium.Icon(color='blue', icon='info_sign')
).add_to(m) # Marker for origin point

folium.Marker(
    location=destination_coords, 
    icon=folium.Icon(color='green', icon='info_sign')
).add_to(m) # Marker for destination point

route_coords = [(G.nodes[node]['y'], G.nodes[node]['x']) for node in shortest_path]
folium.PolyLine(route_coords, color='red', weight=4, opacity=0.8).add_to(m) # Add the route line

# Display the map
# ==============================================================================
m
 Shortest Path: 17558.97 methers
Out[4]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Closest Facility AnalysisΒΆ

In this section, we simulate delivery locations and assign each one to its nearest logistics hub based on road network distance. This is an example of a many-to-one routing problem, commonly used in last-mile delivery planning.

Key Objectives:

  • Generate delivery points within the geographic boundary of logistics centers.
  • Use nearest neighbor lookup and Dijkstra’s algorithm to determine the optimal warehouse for each point.
  • Map the results interactively to visualize how logistics centers serve different delivery zones.

This technique is crucial for logistics companies aiming to balance delivery workloads and reduce response times by choosing the most accessible fulfillment center.

InΒ [5]:
from shapely.geometry import Point

# Get bounding box coordinates of all logistics centers in Rotterdam
# (used to generate random delivery points within the city limits)
# ==============================================================================
minx, miny, maxx, maxy = rotterdam.total_bounds

# Generate 10 random points within the bounding box
# ==============================================================================
num_points = 10
x_coords = np.random.uniform(minx, maxx, num_points)
y_coords = np.random.uniform(miny, maxy, num_points)
delivery_points = [Point(x, y) for x, y in zip(x_coords, y_coords)]

# Create a GeoDataFrame for delivery points
# ==============================================================================
deliveries = gpd.GeoDataFrame(geometry=delivery_points, crs=rotterdam.crs)

# Initialize the map centered around Rotterdam
# ==============================================================================
m = folium.Map(location=[(miny + maxy) / 2, (minx + maxx) / 2], zoom_start=12)

# Build a KDTree of logistics center node coordinates
# ==============================================================================
center_nodes = rotterdam['nearest_node'].values
center_coords = np.array([(G.nodes[node]['x'], G.nodes[node]['y']) for node in center_nodes])
center_tree = cKDTree(center_coords)

# Assign nearest road node to each delivery point
# ==============================================================================
delivery_coords = np.array([(point.x, point.y) for point in deliveries.geometry])
node_ids = list(G.nodes)
node_coords = np.array([(G.nodes[node]['x'], G.nodes[node]['y']) for node in node_ids])
node_tree = cKDTree(node_coords)

# Match each delivery point to its nearest road node
# ==============================================================================
_, node_indices = node_tree.query(delivery_coords, k=1)
deliveries['delivery_node'] = [node_ids[i] for i in node_indices]

# Match each delivery to the closest logistics center based on network distance
# ==============================================================================
delivery_routes = []
for delivery_node in deliveries['delivery_node']:
    # Use Dijkstra to all centers
    lengths = nx.single_source_dijkstra_path_length(G, delivery_node, weight='length')
    
    # Pick closest center
    min_dist = float('inf')
    closest_center = None
    for center_node in center_nodes:
        if center_node in lengths and lengths[center_node] < min_dist:
            min_dist = lengths[center_node]
            closest_center = center_node
    
    delivery_routes.append({
        'delivery_node': delivery_node,
        'center_node': closest_center,
        'distance': min_dist
    })

# Add logistics center to the map
# ==============================================================================
for idx, row in rotterdam.iterrows():
    center_coords = (row.geometry.centroid.y, row.geometry.centroid.x)
    folium.Marker(location=center_coords, icon=folium.Icon(color='blue', icon='industry')).add_to(m)

# Add delivery points and routes to the map
# ==============================================================================
for idx, row in deliveries.iterrows():
    delivery_coords = (row.geometry.y, row.geometry.x)
    folium.Marker(location=delivery_coords, icon=folium.Icon(color='green', icon='home')).add_to(m)
    
    # Find the corresponding route
    route_info = delivery_routes[idx]
    path = nx.shortest_path(G, route_info['delivery_node'], route_info['center_node'], weight='length')
    route_coords = [(G.nodes[node]['y'], G.nodes[node]['x']) for node in path]
    folium.PolyLine(route_coords, color='red', weight=2.5, opacity=0.8).add_to(m)


# Display the map
# ==============================================================================
m
Out[5]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Service Area AnalysisΒΆ

The Service Area Analysis estimates the geographic coverage or reachability from a logistics center within a certain distance threshold (e.g., 10 km of travel distance). This technique is essential in:

  • Territory planning (e.g., determining coverage areas of warehouses)
  • Accessibility studies (e.g., identifying underserved areas)
  • Network resilience planning (e.g., evaluating which roads are critical to maintain access)
InΒ [9]:
# Choose a few random logistics centers (Origin Nodes)
# ==============================================================================
origin_nodes = rotterdam.iloc[[5, 24, 70]]['nearest_node'].values

colors = ['red', 'green', 'blue']
service_dists = [10000]  # or multiple if you want [5000, 10000]

m = folium.Map(location=origin_coords, zoom_start=11)

# Loop over each logistics center to generate and visualize its 10 km service area
# ==================================================================================
for idx, origin_node in enumerate(origin_nodes):
    for dist in service_dists:
        subgraph = nx.ego_graph(G, origin_node, radius=dist, distance='length')
        edges = ox.graph_to_gdfs(subgraph, nodes=False, edges=True)
        
        for _, row in edges.iterrows():
            line = folium.PolyLine(
                locations=[(y, x) for x, y in row['geometry'].coords],
                color=colors[idx],
                weight=2.5,
                opacity=0.5
            )
            line.add_to(m)

    # Add a marker for the center
    # ==============================================================================
    origin_coords = (G.nodes[origin_node]['y'], G.nodes[origin_node]['x'])
    folium.Marker(location=origin_coords, icon=folium.Icon(color=colors[idx], icon='industry')).add_to(m)

m
Out[9]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Facility Location OptimizationΒΆ

In this section, we tackled the classic Facility Location Optimization (FLO) problem from a logistics and supply chain lens: Where should we place warehouses (facilities) to best serve a set of delivery points (customers)?

Our aim was to:

  • Open a limited number of facilities from a candidate pool.
  • Assign each delivery point to one open facility.
  • Minimize the total distance from all delivery points to their assigned facilities.

Although we previously worked with a GeoDataFrame named rotterdam containing actual logistics buildings and attributes (e.g., footprint, job density), we intentionally avoided using it here for three key reasons:

  1. Simplification for Learning: The rotterdam data included a lot of extra variables and spatial geometry. For our first implementation of FLO, we focused on a cleaner, minimal dataset to make the optimization logic easier to follow.
  2. Flexibility & Randomization: Using graph nodes directly gave us flexibility to randomly generate different scenarios and test the optimization process under various configurations.
  3. Avoiding CRS/Geometry Issues: The rotterdam dataframe has polygon geometries, whereas OSMnx graph nodes are points. Mixing them would require converting polygons to centroids and ensuring coordinate reference systems (CRS) matched, which introduces complexity. For clarity and faster execution, we kept everything at the node level.
InΒ [20]:
import osmnx as ox
from pulp import *

# Download the Road Network
# ==============================================================================
city = "Rotterdam, Netherlands"
G = ox.graph_from_place(city, network_type="drive")
nodes, edges = ox.graph_to_gdfs(G)

np.random.seed(42)

# Sample demand points (e.g., customers)
# ==============================================================================
demand_points = nodes.sample(10).reset_index()
demand_points['id'] = ['D' + str(i) for i in range(len(demand_points))]

# Sample facility candidates (e.g., warehouses)
# ==============================================================================
facility_points = nodes.sample(5).reset_index()
facility_points['id'] = ['F' + str(i) for i in range(len(facility_points))]

# Calculate Distance Matrix
# ==============================================================================
demand_coords = demand_points[['y', 'x']].to_numpy()
facility_coords = facility_points[['y', 'x']].to_numpy()

distance_matrix = np.round(
    np.linalg.norm(demand_coords[:, None] - facility_coords, axis=2), 2
)

# Build and Solve the Optimization Model
# ==============================================================================
model = LpProblem("Facility_Location_Problem", LpMinimize)

# Decision variable: assign[d][f] = 1 if demand d assigned to facility f
# ==============================================================================
assign = LpVariable.dicts("Assign", 
                          [(d, f) for d in range(len(demand_points)) for f in range(len(facility_points))],
                          cat='Binary')

# Open facility variable: 1 if facility f is open
# ==============================================================================
open_facility = LpVariable.dicts("Open", 
                                 list(range(len(facility_points))), 
                                 cat='Binary')

# Objective: Minimize total distance
# ==============================================================================
model += lpSum(assign[d, f] * distance_matrix[d][f]
               for d in range(len(demand_points))
               for f in range(len(facility_points)))

# Constraint: Each demand must be assigned to one facility
# ==============================================================================
for d in range(len(demand_points)):
    model += lpSum(assign[d, f] for f in range(len(facility_points))) == 1

# Constraint: Only assign to open facilities
# ==============================================================================
for d in range(len(demand_points)):
    for f in range(len(facility_points)):
        model += assign[d, f] <= open_facility[f]

# Limit number of open facilities
# ==============================================================================
model += lpSum(open_facility[f] for f in range(len(facility_points))) <= 2

# Solve it
# ==============================================================================
model.solve()
Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /opt/anaconda3/lib/python3.12/site-packages/pulp/apis/../solverdir/cbc/osx/i64/cbc /var/folders/p4/km2q7qsd1cl5b1h60220jf8w0000gn/T/f29e0b4fc3e04cc7a473591e24688089-pulp.mps -timeMode elapsed -branch -printingOptions all -solution /var/folders/p4/km2q7qsd1cl5b1h60220jf8w0000gn/T/f29e0b4fc3e04cc7a473591e24688089-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 66 COLUMNS
At line 382 RHS
At line 444 BOUNDS
At line 500 ENDATA
Problem MODEL has 61 rows, 55 columns and 155 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 0.36 - 0.00 seconds
Cgl0004I processed model has 61 rows, 55 columns (55 integer (55 of which binary)) and 155 elements
Cutoff increment increased from 1e-05 to 0.00999
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 0.36
Cbc0038I Before mini branch and bound, 55 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (0.01 seconds)
Cbc0038I After 0.01 seconds - Feasibility pump exiting with objective of 0.36 - took 0.00 seconds
Cbc0012I Integer solution of 0.36 found by feasibility pump after 0 iterations and 0 nodes (0.01 seconds)
Cbc0001I Search completed - best objective 0.36, took 0 iterations and 0 nodes (0.01 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 0.36 to 0.36
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                0.36000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.02

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.02

Out[20]:
1
InΒ [22]:
# Extract coordinates
# ==============================================================================
demand_coords = demand_points[['y', 'x']].to_numpy()
facility_coords = facility_points[['y', 'x']].to_numpy()

# Identify which facilities are open
# ==============================================================================
open_facilities = [f for f in range(len(facility_points)) if open_facility[f].varValue == 1]

# Create a plot
# ==============================================================================
fig, ax = plt.subplots(figsize=(10, 10))

# Plot the road network as background
# ==============================================================================
edges.plot(ax=ax, linewidth=0.2, color='gray', alpha=0.4)

# Plot demand points
# ==============================================================================
plt.scatter(demand_coords[:, 1], demand_coords[:, 0], c='red', label='Demand Points')

# Plot opened facilities
# ==============================================================================
for f in open_facilities:
    plt.scatter(facility_coords[f][1], facility_coords[f][0], c='blue', s=100, label='Opened Facility' if f == open_facilities[0] else "")

# Draw assignment lines
# ==============================================================================
for d in range(len(demand_points)):
    for f in open_facilities:
        if assign[d, f].varValue == 1:
            x = [demand_coords[d][1], facility_coords[f][1]]
            y = [demand_coords[d][0], facility_coords[f][0]]
            plt.plot(x, y, 'g--', linewidth=1)

# Add legend and title
# ==============================================================================
plt.legend()
plt.title("Facility Location Optimization: Opened Facilities and Assignments")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.grid(True)
plt.show()
No description has been provided for this image

ConclusionsΒΆ

  1. Real-World Network Data Enables Richer Logistics Insights

By leveraging real road network data from OSMnx, we moved beyond theoretical models and worked with geographically accurate spatial relationships. This approach allowed us to:

  • Compute shortest paths based on real travel distances.
  • Understand service areas using realistic reachability constraints.
  • Identify optimal facility locations in a way that reflects actual urban infrastructure.
  1. Optimization Provides Data-Driven Warehouse Siting

Through Facility Location Optimization, we demonstrated how a logistics company might:

  • Minimize overall travel distance or cost,
  • Choose a limited number of warehouses to open,
  • Assign delivery locations efficiently based on proximity.

Even with a simplified model (Euclidean distances, no capacity constraints), this approach generated an actionable, interpretable assignment plan that reduced delivery inefficiencies.

  1. Service Area & Shortest Route Analysis Add Operational Value

We visualized:

  • Shortest paths between selected warehouses,
  • Service areas within 5 km and 10 km radii,
  • Closest facility assignments for random delivery points.

These analyses are not only useful for strategic planning (where to build), but also for daily operations (which warehouse should serve which customer).

  1. Geospatial Optimization Bridges the Gap Between Mapping and Modeling

This project showed how geospatial analysis (using tools like GeoPandas, OSMnx, and Folium) can be integrated with mathematical optimization (via PuLP) to solve complex urban logistics challenges in a replicable, scalable way.

  1. Next Steps and Real-World Applications

To move this from prototype to production, one could:

  • Incorporate real customer delivery data and existing warehouse locations (like those in the rotterdam GeoDataFrame),
  • Replace Euclidean distances with actual road travel times,
  • Add constraints like facility capacities, delivery time windows, or environmental impact.